Drawing vector maps with ggplot2 and sf

library(tidyverse)
library(sf)

options(digits = 3)
set.seed(1234)
theme_set(theme_minimal())

Unlike raster image maps, vector maps require you to obtain spatial data files which contain detailed information necessary to draw all the components of a map (e.g. points, lines, polygons). Once you successfully import that data into R, ggplot2 works with simple features data frames to easily generate geospatial visualizations using all the core elements and approaches of ggplot().

Import USA state boundaries

First we will import a spatial data file containing the boundaries of all 50 states in the United States¹ using sf::st_read():

usa <- st_read("data/census_bureau/cb_2013_us_state_20m/cb_2013_us_state_20m.shp")

## Reading layer `cb_2013_us_state_20m' from data source `/Users/soltoffbc/Projects/Computing for Social Sciences/uc-cfss.github.io/data/census_bureau/cb_2013_us_state_20m/cb_2013_us_state_20m.shp' using driver `ESRI Shapefile'
## Simple feature collection with 52 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179 ymin: 17.9 xmax: 180 ymax: 71.4
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs

Draw the boundaries

ggplot2 contains a geometric object specifically for simple feature objects called geom_sf(). This works reasonably well when you need to draw polygons, like our state boundaries. Support for simple features in ggplot2 is under active development, so you may not find adequate support for plotting line or point features. To draw the map, we pass the simple features data frame as the data argument.

ggplot(data = usa) +
  geom_sf()

Because simple features data frames are standardized with the geometry column always containing information on the geographic coordinates of the features, we do not need to specify additional parameters for aes(). Notice a problem with the map above: it wastes a lot of space. This is caused by the presence of Alaska and Hawaii in the dataset. The Aleutian Islands cross the the 180th meridian, requiring the map to show the Eastern hemisphere. Likewise, Hawaii is substantially distant from the continental United States.

Plot a subset of a map

One solution is to plot just the lower 48 states. That is, exclude Alaska and Hawaii, as well as DC and Puerto Rico.² Because simple features data frames contain one row per feature and in this example a feature is defined as a state, we can use filter() from dplyr to exclude these four states/territories.

(usa_48 <- usa %>%
  filter(!(NAME %in% c("Alaska", "District of Columbia", "Hawaii", "Puerto Rico"))))

## Simple feature collection with 48 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -125 ymin: 24.5 xmax: -66.9 ymax: 49.4
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
## First 10 features:
##    STATEFP  STATENS    AFFGEOID GEOID STUSPS        NAME LSAD    ALAND
## 1       01 01779775 0400000US01    01     AL     Alabama   00 1.31e+11
## 2       05 00068085 0400000US05    05     AR    Arkansas   00 1.35e+11
## 3       06 01779778 0400000US06    06     CA  California   00 4.03e+11
## 4       09 01779780 0400000US09    09     CT Connecticut   00 1.25e+10
## 5       12 00294478 0400000US12    12     FL     Florida   00 1.39e+11
## 6       13 01705317 0400000US13    13     GA     Georgia   00 1.49e+11
## 7       16 01779783 0400000US16    16     ID       Idaho   00 2.14e+11
## 8       17 01779784 0400000US17    17     IL    Illinois   00 1.44e+11
## 9       18 00448508 0400000US18    18     IN     Indiana   00 9.28e+10
## 10      20 00481813 0400000US20    20     KS      Kansas   00 2.12e+11
##      AWATER                       geometry
## 1  4.59e+09 MULTIPOLYGON (((-88.3 30.2,...
## 2  2.96e+09 MULTIPOLYGON (((-94.6 36.5,...
## 3  2.05e+10 MULTIPOLYGON (((-119 33.5, ...
## 4  1.82e+09 MULTIPOLYGON (((-73.7 41.1,...
## 5  3.14e+10 MULTIPOLYGON (((-80.7 24.9,...
## 6  4.95e+09 MULTIPOLYGON (((-85.6 35, -...
## 7  2.40e+09 MULTIPOLYGON (((-117 44.4, ...
## 8  6.20e+09 MULTIPOLYGON (((-91.5 40.2,...
## 9  1.54e+09 MULTIPOLYGON (((-88.1 37.9,...
## 10 1.35e+09 MULTIPOLYGON (((-102 40, -1...

ggplot(data = usa_48) +
  geom_sf()

Since the map is a ggplot() object, it can easily be modified like any other ggplot() graph. We could change the color of the map and the borders:

ggplot(data = usa_48) +
  geom_sf(fill = "palegreen", color = "black")

fiftystater

Rather than excluding them entirely, most maps of the United States place Alaska and Hawaii as insets to the south of California. Until recently, in R this was an extremely tedious task that required manually changing the latitude and longitude coordinates for these states to place them in the correct location. Fortunately a new package is available that has already done the work for you. fiftystater includes the fifty_states data frame which contains adjusted coordinates for Alaska and Hawaii to plot them with the mainland. It can be installed from GitHub using devtools::install_github("wmurphyrd/fiftystater").

fifty_states is an ordinary data frame with one row for each point:

library(fiftystater)

data("fifty_states")
as_tibble(fifty_states)

## # A tibble: 13,694 x 7
##     long   lat order hole  piece id      group    
##    <dbl> <dbl> <int> <lgl> <fct> <chr>   <fct>    
##  1 -85.1  32.0     1 FALSE 1     alabama Alabama.1
##  2 -85.1  31.9     2 FALSE 1     alabama Alabama.1
##  3 -85.1  31.9     3 FALSE 1     alabama Alabama.1
##  4 -85.1  31.8     4 FALSE 1     alabama Alabama.1
##  5 -85.1  31.8     5 FALSE 1     alabama Alabama.1
##  6 -85.1  31.7     6 FALSE 1     alabama Alabama.1
##  7 -85.1  31.7     7 FALSE 1     alabama Alabama.1
##  8 -85.1  31.7     8 FALSE 1     alabama Alabama.1
##  9 -85.1  31.6     9 FALSE 1     alabama Alabama.1
## 10 -85.0  31.6    10 FALSE 1     alabama Alabama.1
## # ... with 13,684 more rows

That is, each row defines a boundary point. The boundary for a state is created by drawing a line connecting each point. We need to convert this points-based data frame to a simple features data frame where each feature is a state.³ First we convert the data frame to a simple features data frame using st_as_sf(). This creates a simple features data frame where each feature is a boundary point, like before.

st_as_sf(fifty_states, coords = c("long", "lat"))

## Simple feature collection with 13694 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -125 ymin: 23.5 xmax: -67 ymax: 49.4
## epsg (SRID):    NA
## proj4string:    NA
## First 10 features:
##    order  hole piece      id     group           geometry
## 1      1 FALSE     1 alabama Alabama.1   POINT (-85.1 32)
## 2      2 FALSE     1 alabama Alabama.1 POINT (-85.1 31.9)
## 3      3 FALSE     1 alabama Alabama.1 POINT (-85.1 31.9)
## 4      4 FALSE     1 alabama Alabama.1 POINT (-85.1 31.8)
## 5      5 FALSE     1 alabama Alabama.1 POINT (-85.1 31.8)
## 6      6 FALSE     1 alabama Alabama.1 POINT (-85.1 31.7)
## 7      7 FALSE     1 alabama Alabama.1 POINT (-85.1 31.7)
## 8      8 FALSE     1 alabama Alabama.1 POINT (-85.1 31.7)
## 9      9 FALSE     1 alabama Alabama.1 POINT (-85.1 31.6)
## 10    10 FALSE     1 alabama Alabama.1   POINT (-85 31.6)

Next we group the simple features data frame by id and piece (these variables define each distinct land mass within a state), combine them together using summarize(), and convert them to polygons:

st_as_sf(fifty_states, coords = c("long", "lat")) %>% 
  # convert sets of points to polygons
  group_by(id, piece) %>% 
  summarize(do_union = FALSE) %>%
  st_cast("POLYGON")

## Simple feature collection with 134 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -125 ymin: 23.5 xmax: -67 ymax: 49.4
## epsg (SRID):    NA
## proj4string:    NA
## First 10 features:
##         id piece do_union                       geometry
## 1  alabama     1    FALSE POLYGON ((-85.1 32, -85.1 3...
## 2   alaska     1    FALSE POLYGON ((-120 25.8, -120 2...
## 3   alaska     2    FALSE POLYGON ((-118 25.2, -118 2...
## 4   alaska     3    FALSE POLYGON ((-113 25.4, -113 2...
## 5   alaska     4    FALSE POLYGON ((-114 26.1, -114 2...
## 6   alaska     5    FALSE POLYGON ((-121 27.9, -121 2...
## 7   alaska     6    FALSE POLYGON ((-114 26, -114 26,...
## 8   alaska     7    FALSE POLYGON ((-120 26.6, -120 2...
## 9   alaska     8    FALSE POLYGON ((-120 24.2, -120 2...
## 10  alaska     9    FALSE POLYGON ((-114 25.5, -114 2...

Finally, we combine together all the pieces of land comprising each state. For example, Alabama is a single contiguous landmass, whereas Alaska is 32 distinct pieces of land. This is another group_by() %>% summarize() operation, but only aggregating on id.

# convert fifty_states to an sf data frame
(sf_fifty <- st_as_sf(fifty_states, coords = c("long", "lat")) %>% 
   # convert sets of points to polygons
   group_by(id, piece) %>% 
   summarize(do_union = FALSE) %>%
   st_cast("POLYGON") %>%
   # convert polygons to multipolygons for states with discontinuous regions
   group_by(id) %>%
   summarize())

## Simple feature collection with 51 features and 1 field
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: -125 ymin: 23.5 xmax: -67 ymax: 49.4
## epsg (SRID):    NA
## proj4string:    NA
## # A tibble: 51 x 2
##    id                                                             geometry
##    <chr>                                                        <GEOMETRY>
##  1 alabama              POLYGON ((-85.1 32, -85.1 31.9, -85.1 31.9, -85.1…
##  2 alaska               MULTIPOLYGON (((-123 23.6, -123 23.6, -123 23.7, …
##  3 arizona              POLYGON ((-115 33, -115 33, -115 33, -115 33, -11…
##  4 arkansas             POLYGON ((-94.5 34.2, -94.5 34.5, -94.4 34.7, -94…
##  5 california           MULTIPOLYGON (((-118 32.8, -118 32.8, -119 32.9, …
##  6 colorado             POLYGON ((-102 37.6, -102 37.4, -102 37, -103 37,…
##  7 connecticut          POLYGON ((-73.5 41.5, -73.5 41.7, -73.5 42, -73 4…
##  8 delaware             POLYGON ((-75.7 38.6, -75.7 38.6, -75.7 38.8, -75…
##  9 district of columbia POLYGON ((-77 39, -76.9 38.9, -77 38.8, -77 38.8,…
## 10 florida              MULTIPOLYGON (((-80.2 25.4, -80.4 25.2, -80.6 25,…
## # ... with 41 more rows

Finally we set the projection method⁴ and draw the map:

st_crs(sf_fifty) <- 4326

ggplot(data = sf_fifty) +
  geom_sf()

Add data to the map

Region boundaries serve as the background in geospatial data visualization - so now we need to add data. Some types of geographic data (points and symbols) are overlaid on top of the boundaries, whereas other data (fill) are incorporated into the region layer itself.

Points

Let’s use our usa_48 map data to add some points. The airports data frame in the nycflights13 package includes geographic info on airports in the United States.

library(nycflights13)
airports

## # A tibble: 1,458 x 8
##    faa   name                   lat    lon   alt    tz dst   tzone        
##    <chr> <chr>                <dbl>  <dbl> <int> <dbl> <chr> <chr>        
##  1 04G   Lansdowne Airport     41.1  -80.6  1044    -5 A     America/New_…
##  2 06A   Moton Field Municip…  32.5  -85.7   264    -6 A     America/Chic…
##  3 06C   Schaumburg Regional   42.0  -88.1   801    -6 A     America/Chic…
##  4 06N   Randall Airport       41.4  -74.4   523    -5 A     America/New_…
##  5 09J   Jekyll Island Airpo…  31.1  -81.4    11    -5 A     America/New_…
##  6 0A9   Elizabethton Munici…  36.4  -82.2  1593    -5 A     America/New_…
##  7 0G6   Williams County Air…  41.5  -84.5   730    -5 A     America/New_…
##  8 0G7   Finger Lakes Region…  42.9  -76.8   492    -5 A     America/New_…
##  9 0P2   Shoestring Aviation…  39.8  -76.6  1000    -5 U     America/New_…
## 10 0S9   Jefferson County In…  48.1 -123.    108    -8 A     America/Los_…
## # ... with 1,448 more rows

Each airport has it’s geographic location encoded through lat and lon. To draw these points on the map, basically we draw a scatterplot with x = lon and y = lat. In fact we could simply do that:

ggplot(airports, aes(lon, lat)) +
  geom_point()

Let’s overlay it with the mapped state borders:

ggplot(data = usa_48) + 
  geom_sf() + 
  geom_point(data = airports, aes(x = lon, y = lat), shape = 1)

Slight problem. We have airports listed outside of the continental United States. There are a couple ways to rectify this. Unfortunately airports does not include a variable identifying state so the filter() operation is not that simple. The easiest solution is to crop the limits of the graph using coord_sf() to only show the mainland:

ggplot(data = usa_48) + 
  geom_sf() + 
  geom_point(data = airports, aes(x = lon, y = lat), shape = 1) +
  coord_sf(xlim = c(-130, -60),
           ylim = c(20, 50))

Alternatively, we can use st_as_sf() to convert airports to a simple features data frame.

airports_sf <- st_as_sf(airports, coords = c("lon", "lat"))
st_crs(airports_sf) <- 4326   # set the coordinate reference system
airports_sf

## Simple feature collection with 1458 features and 6 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -177 ymin: 19.7 xmax: 174 ymax: 72.3
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 1,458 x 7
##    faa   name              alt    tz dst   tzone                  geometry
##    <chr> <chr>           <int> <dbl> <chr> <chr>               <POINT [°]>
##  1 04G   Lansdowne Airp…  1044    -5 A     America/N…         (-80.6 41.1)
##  2 06A   Moton Field Mu…   264    -6 A     America/C…         (-85.7 32.5)
##  3 06C   Schaumburg Reg…   801    -6 A     America/C…           (-88.1 42)
##  4 06N   Randall Airport   523    -5 A     America/N…         (-74.4 41.4)
##  5 09J   Jekyll Island …    11    -5 A     America/N…         (-81.4 31.1)
##  6 0A9   Elizabethton M…  1593    -5 A     America/N…         (-82.2 36.4)
##  7 0G6   Williams Count…   730    -5 A     America/N…         (-84.5 41.5)
##  8 0G7   Finger Lakes R…   492    -5 A     America/N…         (-76.8 42.9)
##  9 0P2   Shoestring Avi…  1000    -5 U     America/N…         (-76.6 39.8)
## 10 0S9   Jefferson Coun…   108    -8 A     America/L…          (-123 48.1)
## # ... with 1,448 more rows

coords tells st_as_sf() which columns contain the geographic coordinates of each airport. To graph the points on the map, we use a second geom_sf():

ggplot() + 
  geom_sf(data = usa_48) + 
  geom_sf(data = airports_sf, shape = 1) +
  coord_sf(xlim = c(-130, -60),
           ylim = c(20, 50))

Symbols

We can change the size or type of symbols on the map. For instance, we can draw a bubble plot (also known as a proportional symbol map) and encode the altitude of the airport through the size channel:

ggplot(data = usa_48) + 
  geom_sf() + 
  geom_point(data = airports, aes(x = lon, y = lat, size = alt),
             fill = "grey", color = "black", alpha = .2) +
  coord_sf(xlim = c(-130, -60),
           ylim = c(20, 50)) +
  scale_size_area(guide = FALSE)

Circle area is proportional to the airport’s altitude (in feet). Or we could scale it based on the number of arriving flights in flights:

airports_n <- flights %>%
  count(dest) %>%
  left_join(airports, by = c("dest" = "faa"))

ggplot(data = usa_48) + 
  geom_sf() + 
  geom_point(data = airports_n, aes(x = lon, y = lat, size = n),
             fill = "grey", color = "black", alpha = .2) +
  coord_sf(xlim = c(-130, -60),
           ylim = c(20, 50)) +
  scale_size_area(guide = FALSE)

airports contains a list of virtually all commercial airports in the United States. However flights only contains data on flights departing from New York City airports (JFK, LaGuardia, or Newark) and only services a few airports around the country.

Fill (choropleths)

Choropleth maps encode information by assigning shades of colors to defined areas on a map (e.g. countries, states, counties, zip codes). There are lots of ways to tweak and customize these graphs, which is generally a good idea because remember that color is one of the harder-to-decode channels.

We will continue to use the usa_48 simple features data frame and draw a choropleth for the number of foreign-born individuals in each state. We get those files from the census_bureau folder. Let’s also normalize our measure by the total population to get the rate of foreign-born individuals in the population:

(fb_state <- read_csv("data/census_bureau/ACS_13_5YR_B05012_state/ACS_13_5YR_B05012.csv") %>%
  mutate(rate = HD01_VD03 / HD01_VD01))

## # A tibble: 51 x 10
##    GEO.id GEO.id2 `GEO.display-la… HD01_VD01 HD02_VD01 HD01_VD02 HD02_VD02
##    <chr>  <chr>   <chr>                <int> <chr>         <int>     <int>
##  1 04000… 01      Alabama            4799277 <NA>        4631045      2881
##  2 04000… 02      Alaska              720316 <NA>         669941      1262
##  3 04000… 04      Arizona            6479703 <NA>        5609835      7725
##  4 04000… 05      Arkansas           2933369 <NA>        2799972      2568
##  5 04000… 06      California        37659181 <NA>       27483342     30666
##  6 04000… 08      Colorado           5119329 <NA>        4623809      5778
##  7 04000… 09      Connecticut        3583561 <NA>        3096374      5553
##  8 04000… 10      Delaware            908446 <NA>         831683      2039
##  9 04000… 11      District of Col…    619371 <NA>         534142      2017
## 10 04000… 12      Florida           19091156 <NA>       15392410     16848
## # ... with 41 more rows, and 3 more variables: HD01_VD03 <int>,
## #   HD02_VD03 <int>, rate <dbl>

Join the data

Now that we have our data, we want to draw it on the map. fb_state contains one row per state, as does usa_48. Since there is a one-to-one match between the data frames, we join the data frames together first, then use that single data frame to draw the map. This differs from the approach above for drawing points because a point feature is not the same thing as a polygon feature. That is, there were more airports then there were states. Because the spatial data is stored in a data frame with one row per state, all we need to do is merge the data frames together on a column that uniquely identifies each row in each data frame.

(usa_fb <- usa_48 %>%
  left_join(fb_state, by = c("STATEFP" = "GEO.id2")))

## Simple feature collection with 48 features and 18 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -125 ymin: 24.5 xmax: -66.9 ymax: 49.4
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
## First 10 features:
##    STATEFP  STATENS    AFFGEOID GEOID STUSPS        NAME LSAD    ALAND
## 1       01 01779775 0400000US01    01     AL     Alabama   00 1.31e+11
## 2       05 00068085 0400000US05    05     AR    Arkansas   00 1.35e+11
## 3       06 01779778 0400000US06    06     CA  California   00 4.03e+11
## 4       09 01779780 0400000US09    09     CT Connecticut   00 1.25e+10
## 5       12 00294478 0400000US12    12     FL     Florida   00 1.39e+11
## 6       13 01705317 0400000US13    13     GA     Georgia   00 1.49e+11
## 7       16 01779783 0400000US16    16     ID       Idaho   00 2.14e+11
## 8       17 01779784 0400000US17    17     IL    Illinois   00 1.44e+11
## 9       18 00448508 0400000US18    18     IN     Indiana   00 9.28e+10
## 10      20 00481813 0400000US20    20     KS      Kansas   00 2.12e+11
##      AWATER      GEO.id GEO.display-label HD01_VD01 HD02_VD01 HD01_VD02
## 1  4.59e+09 0400000US01           Alabama   4799277      <NA>   4631045
## 2  2.96e+09 0400000US05          Arkansas   2933369      <NA>   2799972
## 3  2.05e+10 0400000US06        California  37659181      <NA>  27483342
## 4  1.82e+09 0400000US09       Connecticut   3583561      <NA>   3096374
## 5  3.14e+10 0400000US12           Florida  19091156      <NA>  15392410
## 6  4.95e+09 0400000US13           Georgia   9810417      <NA>   8859747
## 7  2.40e+09 0400000US16             Idaho   1583364      <NA>   1489560
## 8  6.20e+09 0400000US17          Illinois  12848554      <NA>  11073828
## 9  1.54e+09 0400000US18           Indiana   6514861      <NA>   6206801
## 10 1.35e+09 0400000US20            Kansas   2868107      <NA>   2677007
##    HD02_VD02 HD01_VD03 HD02_VD03   rate                       geometry
## 1       2881    168232      2881 0.0351 MULTIPOLYGON (((-88.3 30.2,...
## 2       2568    133397      2568 0.0455 MULTIPOLYGON (((-94.6 36.5,...
## 3      30666  10175839     30666 0.2702 MULTIPOLYGON (((-119 33.5, ...
## 4       5553    487187      5553 0.1360 MULTIPOLYGON (((-73.7 41.1,...
## 5      16848   3698746     16848 0.1937 MULTIPOLYGON (((-80.7 24.9,...
## 6       7988    950670      7988 0.0969 MULTIPOLYGON (((-85.6 35, -...
## 7       2528     93804      2528 0.0592 MULTIPOLYGON (((-117 44.4, ...
## 8      10091   1774726     10093 0.1381 MULTIPOLYGON (((-91.5 40.2,...
## 9       4499    308060      4500 0.0473 MULTIPOLYGON (((-88.1 37.9,...
## 10      3095    191100      3100 0.0666 MULTIPOLYGON (((-102 40, -1...

Draw the map

With the newly combined data frame, use geom_sf() and define the fill aesthetic based on the column in usa_fb you want to visualize.

ggplot(data = usa_fb) +
  geom_sf(aes(fill = rate))

Bin data to discrete intervals

When creating a heatmap with a continuous variable, one must decide whether to keep the variable as continuous or collapse it into a series of bins with discrete colors. While keep the variable continuous is technically more precise, the human eye cannot usually distinguish between two colors which are very similar to one another. By converting the variable to a discrete variable, you easily distinguish between the different levels. If you decide to convert a continuous variable to a discrete variable, you will need to decide how to do this. While cut() is a base R function for converting continuous variables into discrete values, ggplot2 offers two functions that explicitly define how we want to bin the numeric vector (column).

cut_interval() makes n groups with equal range:

usa_fb %>%
  mutate(rate_cut = cut_interval(rate, 6)) %>%
  ggplot() +
  geom_sf(aes(fill = rate_cut))

Whereas cut_number() makes n groups with (approximately) equal numbers of observations:

usa_fb %>%
  mutate(rate_cut = cut_number(rate, 6)) %>%
  ggplot() +
  geom_sf(aes(fill = rate_cut))

See this StackOverflow thread for a more in-depth discussion on the merits of bucketizing a continuous variable and whether to use cut_interval() or cut_number().

Changing map projection

Representing portions of the globe on a flat surface can be challenging. Depending on how you project the map, you can distort or emphasize certain features of the map. Fortunately, ggplot() includes the coord_sf() function which allows us to easily implement different projection methods. In order to implement coordinate transformations, you need to know the coordinate reference system that defines the projection method. The “easiest” approach is to provide what is known as the proj4string that defines the projection method. PROJ4 is a generic coordinate transformation software that allows you to convert between projection methods. If you get really into geospatial analysis and visualization, it is helpful to learn this system.

For our purposes here, proj4string is a character string in R that defines the coordinate system and includes parameters specific to a given coordinate transformation. PROJ4 includes some documentation on common projection methods that can get you started. Some projection methods are relatively simple and require just the name of the projection, like for a Mercator projection ("+proj=merc"):

map_proj_base <- ggplot(data = usa_48) +
  geom_sf()

map_proj_base +
  coord_sf(crs = "+proj=merc") +
  ggtitle("Mercator projection")

Other coordinate systems require specification of the standard lines, or lines that define areas of the surface of the map that are tangent to the globe. These include Gall-Peters, Albers equal-area, and Lambert azimuthal.

map_proj_base +
  coord_sf(crs = "+proj=cea +lon_0=0 +lat_ts=45") +
  ggtitle("Gall-Peters projection")

map_proj_base +
  coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100") +
  ggtitle("Albers equal-area projection")

map_proj_base +
  coord_sf(crs = "+proj=laea +lat_0=35 +lon_0=-100") +
  ggtitle("Lambert azimuthal projection")

Session Info

devtools::session_info()

##  setting  value                       
##  version  R version 3.5.1 (2018-07-02)
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/Chicago             
##  date     2019-01-02                  
## 
##  package    * version date       source        
##  assertthat   0.2.0   2017-04-11 CRAN (R 3.5.0)
##  backports    1.1.2   2017-12-13 CRAN (R 3.5.0)
##  base       * 3.5.1   2018-07-05 local         
##  base64enc    0.1-3   2015-07-28 CRAN (R 3.5.0)
##  bindr        0.1.1   2018-03-13 CRAN (R 3.5.0)
##  bindrcpp   * 0.2.2   2018-03-29 CRAN (R 3.5.0)
##  broom        0.5.0   2018-07-17 CRAN (R 3.5.0)
##  cellranger   1.1.0   2016-07-27 CRAN (R 3.5.0)
##  class        7.3-14  2015-08-30 CRAN (R 3.5.1)
##  classInt     0.2-3   2018-04-16 CRAN (R 3.5.0)
##  cli          1.0.0   2017-11-05 CRAN (R 3.5.0)
##  colorspace   1.3-2   2016-12-14 CRAN (R 3.5.0)
##  compiler     3.5.1   2018-07-05 local         
##  crayon       1.3.4   2017-09-16 CRAN (R 3.5.0)
##  datasets   * 3.5.1   2018-07-05 local         
##  DBI          1.0.0   2018-05-02 CRAN (R 3.5.0)
##  devtools     1.13.6  2018-06-27 CRAN (R 3.5.0)
##  digest       0.6.18  2018-10-10 cran (@0.6.18)
##  dplyr      * 0.7.8   2018-11-10 cran (@0.7.8) 
##  e1071        1.7-0   2018-07-28 CRAN (R 3.5.0)
##  evaluate     0.11    2018-07-17 CRAN (R 3.5.0)
##  forcats    * 0.3.0   2018-02-19 CRAN (R 3.5.0)
##  ggplot2    * 3.1.0   2018-10-25 cran (@3.1.0) 
##  glue         1.3.0   2018-07-17 CRAN (R 3.5.0)
##  graphics   * 3.5.1   2018-07-05 local         
##  grDevices  * 3.5.1   2018-07-05 local         
##  grid         3.5.1   2018-07-05 local         
##  gtable       0.2.0   2016-02-26 CRAN (R 3.5.0)
##  haven        1.1.2   2018-06-27 CRAN (R 3.5.0)
##  hms          0.4.2   2018-03-10 CRAN (R 3.5.0)
##  htmltools    0.3.6   2017-04-28 CRAN (R 3.5.0)
##  httr         1.3.1   2017-08-20 CRAN (R 3.5.0)
##  jsonlite     1.5     2017-06-01 CRAN (R 3.5.0)
##  knitr        1.20    2018-02-20 CRAN (R 3.5.0)
##  lattice      0.20-35 2017-03-25 CRAN (R 3.5.1)
##  lazyeval     0.2.1   2017-10-29 CRAN (R 3.5.0)
##  lubridate    1.7.4   2018-04-11 CRAN (R 3.5.0)
##  magrittr     1.5     2014-11-22 CRAN (R 3.5.0)
##  memoise      1.1.0   2017-04-21 CRAN (R 3.5.0)
##  methods    * 3.5.1   2018-07-05 local         
##  modelr       0.1.2   2018-05-11 CRAN (R 3.5.0)
##  munsell      0.5.0   2018-06-12 CRAN (R 3.5.0)
##  nlme         3.1-137 2018-04-07 CRAN (R 3.5.1)
##  pillar       1.3.0   2018-07-14 CRAN (R 3.5.0)
##  pkgconfig    2.0.2   2018-08-16 CRAN (R 3.5.1)
##  plyr         1.8.4   2016-06-08 CRAN (R 3.5.0)
##  purrr      * 0.2.5   2018-05-29 CRAN (R 3.5.0)
##  R6           2.3.0   2018-10-04 cran (@2.3.0) 
##  Rcpp         1.0.0   2018-11-07 cran (@1.0.0) 
##  readr      * 1.1.1   2017-05-16 CRAN (R 3.5.0)
##  readxl       1.1.0   2018-04-20 CRAN (R 3.5.0)
##  rlang        0.3.0.1 2018-10-25 CRAN (R 3.5.0)
##  rmarkdown    1.10    2018-06-11 CRAN (R 3.5.0)
##  rprojroot    1.3-2   2018-01-03 CRAN (R 3.5.0)
##  rstudioapi   0.7     2017-09-07 CRAN (R 3.5.0)
##  rvest        0.3.2   2016-06-17 CRAN (R 3.5.0)
##  scales       1.0.0   2018-08-09 CRAN (R 3.5.0)
##  sf         * 0.7-1   2018-10-24 CRAN (R 3.5.0)
##  spData       0.2.9.3 2018-08-01 CRAN (R 3.5.0)
##  stats      * 3.5.1   2018-07-05 local         
##  stringi      1.2.4   2018-07-20 CRAN (R 3.5.0)
##  stringr    * 1.3.1   2018-05-10 CRAN (R 3.5.0)
##  tibble     * 1.4.2   2018-01-22 CRAN (R 3.5.0)
##  tidyr      * 0.8.1   2018-05-18 CRAN (R 3.5.0)
##  tidyselect   0.2.5   2018-10-11 cran (@0.2.5) 
##  tidyverse  * 1.2.1   2017-11-14 CRAN (R 3.5.0)
##  tools        3.5.1   2018-07-05 local         
##  units        0.6-0   2018-06-09 CRAN (R 3.5.0)
##  utils      * 3.5.1   2018-07-05 local         
##  withr        2.1.2   2018-03-15 CRAN (R 3.5.0)
##  xml2         1.2.0   2018-01-24 CRAN (R 3.5.0)
##  yaml         2.2.0   2018-07-25 CRAN (R 3.5.0)

---

1. Plus the District of Columbia and Puerto Rico↩

2. Issues of political sovereignty aside, these entities are frequently excluded from maps depending on the data to be incorporated. You can always choose to leave them in the map.↩

3. Modified from Convert a fortified data.frame back to sf object on StackOverflow.↩

4. Run ?st_crs from the console to read the documentation for this function and understand how this number is determined.↩

This work is licensed under the CC BY-NC 4.0 Creative Commons License.